#gwa analysis with SNPassoc
###########################

#reading the data
pedigree.ped <- read.table(file="pedigree.dat",header=TRUE,sep=" ")
pedigree.ped[1:10,1:10]
ncol(pedigree.ped)
phenotype.phe <- read.table(file="pheno.dat",header=TRUE,sep=" ")
phenotype.phe[1:10,]

#make subselection of parents only
pedigree.ped <- subset(pedigree.ped, FID==0 & MID==0)
pedigree.ped[1:10,1:10]
nrow(pedigree.ped)
phenotype.phe <- subset(phenotype.phe, phenotype.phe[,1] %in% pedigree.ped[,2])
phenotype.phe[1:10,]
nrow(phenotype.phe)

#pedigree and pheno data preparation for input to SNPassoc functions
library(SNPassoc)

SNPphe <- cbind(pedigree.ped[,7:(ncol(pedigree.ped))],phenotype.phe[,3])
SNPandPHE <- setupSNP(SNPphe, 1:(ncol(SNPphe)-1), sort = FALSE, info, sep = "/") #manual p4
names(SNPandPHE)[1] <- "status"

#perform first quality checks on the data
summary(SNPandPHE)  #manual p5
class(summary(SNPandPHE))
plotMissing(SNPandPHE)  #manual p6
res<-tableHWE(SNPandPHE)   #manual p6
res

#perform a GWA analysis     # manual p9
All<-WGassociation(status~1,data=SNPandPHE,model="all")
summary(All)
print(All)
pvalAll<-pvalues(All)
WGstats(All)
plot(All)
labels(All)

#more elaborate ways to assess genome-wide significance   # manual p12
All.perm<-scanWGassociation(status ~1, data=SNPandPHE,model="log-add", nperm=2)   #take 1000 permutations if possible
res.perm<- permTest(All.perm)
print(res.perm)
plot(res.perm)
res.perm.rtp<- permTest(All.perm,method="rtp",K=20)
print(res.perm.rtp)

#interaction analysis
int.cod <- interactionPval(status~1, data=SNPandPHE,model="codominant")
int.dom <- interactionPval(status~1, data=SNPandPHE,model="dominant")
int.rec <- interactionPval(status~1, data=SNPandPHE,model="recessive")
int.ove <- interactionPval(status~1, data=SNPandPHE,model="overdominant")
int.log <- interactionPval(status~1, data=SNPandPHE,model="log-additive")
print(int.log)
plot(int.log)

#basic random forests approach
library(randomForest)
as.data.frame(SNPandPHE)
set.seed(1969)
data.rf <- randomForest(status ~ ., data=SNPandPHE, mtry=3,importance=TRUE, na.action=na.omit)
print(data.rf)
round(importance(data.rf), 2) # Show "importance" of variables: higher value means more import

#flexible function to perform machine learning applications
library(MLInterfaces)
help(MLearn)
set.seed(1969)
kp = sample(1:nrow(SNPandPHE), size=75*nrow(SNPandPHE)/100)
data.rf = MLearn(status~., data=crabs, randomForestI, kp, ntree=600
impV = getVarImp(data.rf)
plot(impV, n=15)


#the input data consists of families - that is why we restricted to unrelateds (parents) when using functions from SNPassoc
#####################################

# exploration of the family data
################################

#reading the data
pedigree.ped <- read.table(file="pedigree.dat",header=TRUE,sep=" ")
pedigree.ped[1:10,1:10]
ncol(pedigree.ped)
phenotype.phe <- read.table(file="pheno.dat",header=TRUE,sep=" ")
phenotype.phe[1:10,]

library(pedigree)
#help(calcG)

inputG <- data.matrix(pedigree.ped[,7:ncol(pedigree.ped)])
inputG[inputG == "A/A"] <- 0
inputG[inputG == "A/B"] <- 1
inputG[inputG == "B/A"] <- 1
inputG[inputG == "B/B"] <- 2

calcG(inputG)
     
ped <- pedigree.ped[,2:4]
#inbreeding
ord <- orderPed(ped)
ped2 <- ped[order(ord),]
F <- calcInbreeding(ped2)    #gives an error when there are loops
F
#generation numbers
ord <- orderPed(ped)
ped2 <- ped[order(ord),]
(gens <- countGen(ped2))
gens
#nr of offspring
countOff(ped2)
cbind(ped2[1:10,],countOff(ped2)[1:10])


# CAMP data from FBAT
######################

pedigree.ped <- read.table(file="CAMP.ped",header=TRUE,sep=" ",skip=1)
names.ped <- c("PID", "ID", "FID", "MID", as.matrix(read.table(file="CAMP.ped",header=FALSE,sep=" ",nrow=1)))
names(pedigree.ped) <- names.ped
phenotype.phe <- read.table(file="CAMPZ.phe",header=TRUE,sep=" ",skip=1)
names.phe <- c("PID", "ID", as.matrix(read.table(file="CAMPZ.phe",header=FALSE,sep=" ",nrow=1)))
names(phenotype.phe) <- names.phe

library(pedigree)

calcG(inputG)
     
ped <- pedigree.ped[,2:4]
#inbreeding
ord <- orderPed(ped)
ped2 <- ped[order(ord),]
F <- calcInbreeding(ped2)    #gives an error when there are loops
F
#generation numbers
ord <- orderPed(ped)
ped2 <- ped[order(ord),]
(gens <- countGen(ped2))
gens
#nr of offspring
countOff(ped2)
cbind(ped2[1:10,],countOff(ped2)[1:10])
